
library (foreign)
library (survival)
library (arm)
library (pscl)
library(apsrtable)


######################################################
#
#   Negative Binomial for MID Initiation
#   Thursday, July 26 1:55 p.m.
#######################################################




## Laptop   

TermLimitsData<-  read.dta("/Users/Jeff/Dropbox/TermLimits/Replication/Data/TermLimitsAnalysisData.dta")

summary(TermLimitsData)
#TermLimitsData <- subset(TermLimitsData, year >= 1960)
#attach(TermLimitsData)


##  Model   ##

mod1<- glm.nb (countinit ~ lameduck + hawkish + lameduck_hawkish +
                           rivalry + numbord + parliament + gender +
                                   peace + peace2 + peace3,data=TermLimitsData)
summary(mod1)
odTest(mod1)
logLik(mod1)

#Likelihood ratio test of H0: Poisson, as restricted NB model:
#n.b., the distribution of the test-statistic under H0 is non-standard
#e.g., see help(odTest) for details/references#

#Critical value of test statistic at the alpha= 0.05 level: 2.7055 
#Chi-Square Test Statistic =  2.9212 p-value = 0.04371 
#> logLik(mod1)
#'log Lik.' -805.2504 (df=12)




forms <- strsplit(as.character(mod1$call)[2], split="~", fixed=T)
forms <- unlist(lapply(forms, function(x)paste("~", x, sep="")))
count.form <- as.formula(forms[2])

pred.dat <- data.frame(
  lameduck = c(0,0,0,0,0,1,1,1,1,1), 
  hawkish = c(-1,0,1,2,3,-1,0,1,2,3),
  lameduck_hawkish = c(0,0,0,0,0,-1,0,1,2,3),
  rivalry = c(0,0,0,0,0,0,0,0,0,0),
  numbord = c(2,2,2,2,2,2,2,2,2,2),
  parliament = c(0,0,0,0,0,0,0,0,0,0),
  gender = c(0,0,0,0,0,0,0,0,0,0),
  peace = c(5,5,5,5,5,5,5,5,5,5),
  peace2 = c(25,25,25,25,25,25,25,25,25,25),
  peace3 = c(125,125,125,125,125,125,125,125,125,125)
) 




X.count <- model.matrix(count.form, data=pred.dat)

#all.b <- do.call("c", mod1$coef)
set.seed(22)
b.mat <- mvrnorm(1000, mod1$coef, vcov(mod1))

count.preds <- exp(X.count %*% t(b.mat))

count.ci <- t(apply(count.preds, 1, quantile, c(.5,.025,.975)))




count.accountable <- matrix(count.preds[1:5,],ncol=1000)
count.lame <- matrix(count.preds[6:10,],ncol=1000)
count.diff <- count.lame - count.accountable


pred_accountable <- apply(count.accountable,1,mean)
pred_accountable_low <- apply(count.accountable, 1, quantile, probs = c(0.025),  na.rm = TRUE)
pred_accountable_high <- apply(count.accountable, 1, quantile, probs = c(0.975),  na.rm = TRUE)

pred_lame <- apply(count.lame,1,mean)
pred_lame_low <- apply(count.lame, 1, quantile, probs = c(0.025),  na.rm = TRUE)
pred_lame_high <- apply(count.lame, 1, quantile, probs = c(0.975),  na.rm = TRUE)

pred_diff <- apply(count.diff,1,mean)
pred_diff_low <- apply(count.diff, 1, quantile, probs = c(0.025),  na.rm = TRUE)
pred_diff_high <- apply(count.diff, 1, quantile, probs = c(0.975),  na.rm = TRUE)




#########################################
## Generating Graphing Data Set
#########################################

ruler  <-seq(-1,3,by=1)

    
PredictedValues <- data.frame( pred_accountable, pred_accountable_low, pred_accountable_high,
                              pred_lame, pred_lame_low, pred_lame_high,
                              pred_diff, pred_diff_low, pred_diff_high,ruler)

write.dta(PredictedValues,file="/Users/Jeff/Dropbox/TermLimits/Replication/Models/IndexEmpirical/Conditional/MultipleDisputes/PredictedValues.dta")



